//////////////////////////////////////////////////////////////////////////////
//
//  Created by Philip Mulcahy on 02/02/2012.
//  Copyright (c) 2012 Philip Mulcahy. All rights reserved.
//
//  This file is part of the note-recog library.
//
//  note-recog is free software: you can redistribute it and/or modify
//  it under the terms of version 3 of the GNU Lesser General Public License 
//  as published by the Free Software Foundation.
//
//  note-recog is distributed is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with note-recog.  If not, see <http://www.gnu.org/licenses/>.
//
//////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include "Util.h"
#include "FastFourierTransform.h"
#include <cassert>
#include "Spectrum.h"


FastFourierTransform::FastFourierTransform(const std::vector<short>& inputTimeDomain,
                                           float samplingFreq):
  _inOrderFloats(inputTimeDomain.size()), 
  _samplingFreq(samplingFreq),
  _inputTimeDomain(inputTimeDomain),
  _frameCount(inputTimeDomain.size()),
  _log2frameCount(static_cast<size_t>(log2f(_frameCount))),
  _fftSetup(vDSP_create_fftsetup(_log2frameCount, FFT_RADIX2))
{
    assert( _log2frameCount == log2f(_frameCount) );
    assert( 1 << _log2frameCount == _frameCount );
    
    _freqDomainAbs.resize(_frameCount/2);
    
    _cplx.realp = new float[_frameCount/2];
    _cplx.imagp = new float[_frameCount/2];
}


FastFourierTransform::~FastFourierTransform()
{
    vDSP_destroy_fftsetup(_fftSetup);
    delete[] _cplx.realp;
    delete[] _cplx.imagp;
}


Spectrum FastFourierTransform::calculate()
{
    for( size_t i=0; i!=_frameCount; ++i )
    {
        // convert from short to float
        _inOrderFloats[i] = _inputTimeDomain[i];
    }
    
    vDSP_ctoz(
              reinterpret_cast<COMPLEX *>(&(*(_inOrderFloats.begin()))),
              2,
              &_cplx,
              1,
              _frameCount/2);
    
    // Do the FFT
    vDSP_fft_zrip(_fftSetup, &_cplx, 1, _log2frameCount, FFT_FORWARD);
    
    for( size_t i = 0; i != _frameCount; i+=2 )
    {
        const float real = _cplx.realp[i/2];
        const float imag = _cplx.imagp[i/2];
        const float abs = hypotf(real,imag);
        _freqDomainAbs[i/2] = abs;
    }
    
    return Spectrum( this->abs_begin(),
                     this->abs_end(),
                     _samplingFreq );
}


FastFourierTransform::Functor::Functor(float samplingFreq)
: _staging(2048),
  _transform(new FastFourierTransform(_staging, samplingFreq))
{}


FastFourierTransform::Functor::~Functor()
{
    delete _transform;
}


Spectrum FastFourierTransform::Functor::operator()(const std::vector<short>& time_domain) const
{
    _staging = time_domain;
    return _transform->calculate();
}

